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ABSTRACT 

We investigate phenomenological models of star formation and supernova feedback 
in A''-body/SPH simulations of galaxy formation. First, we compare different prescrip- 
tions in the literature for turning cold gas into stars neglecting feedback effects. We 
find that most prescriptions give broadly similar results: the ratio of cold gas to stars 
in the final galaxies is primarily controlled by the range of gas densities where star 
formation is allowed to proceed efficiently. In the absence of feedback, the fraction 
of gas that cools is much too high resulting, for example, in a X-band luminosity 
function that is much brighter than observed. This problem is ameliorated by includ- 
ing a feedback model which either imparts radial kinetic perturbations to galactic gas 
or directly reheats such material and prevents it from cooling for a certain period of 
time. In both these models, a significant fraction of cold gas is heated and expelled 
from haloes with an efficiency that varies inversely with with halo circular velocity. 
Increasing the resolution of a simulation allows a wider dynamic range in mass to be 
followed, but the average properties of the resolved galaxy population remain largely 
unaffected. However, as the resolution is increased, more and more gas is reheated by 
small galaxies; our results suggest that convergence requires the full mass range of 
galaxies to be resolved. 
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1 INTRODUCTION 

The first detailed modeUing; of galaxy fo rmation more than 
twenty years ago ( White fc Rees 1978| ) revealed a signifi- 
cant difficulty in hierarchical clustering cosmologies. These 
assume that galaxies form when gas cools and condenses 
in dark matter haloes. However, at early times, pre-galactic 
gas clouds are so dense that the cooling time is extremely 
short. Consequently, in the absence of other processes, bary- 
onic maUei cools caLasLiuphictilly In sub-galacLic haluys and 



These ideas carry through immediately to modern cos- 
mological the ory, which is based up on the cold dark mat- 
ter paradigm (White & Frenk 1991). This provides a natu- 
ral motivation for hierarchical clustering and so it inherits 
the "cooling catastrophe", generic to this kind of models. 
There is now a growing body of work attempting to un- 
derstand likely s ources of feedback (e.g. tPekcl & Silk 198(: 
Efstathiou 19921 



Thoul fc Weinberg 1996 



Mac Low & Ferrara 199£ : Efstathiou 200C 



Silk & Rees 199J 



Madau, Ferrara 



forms s tars, leaving no gas available to make up the intcr- 



fc Rees 200d|r 

their efi^ec t s on models of galaxy formation 



galactiti medium obsuived aL high ledshifL, ui Lhe inUacIus- 
ter mec 



L ubjLi ULd aL luw itdjIufL. WhiLt &, Rllj piupujud 



White 1993; Kauffmann ct al. 1994|; |Colc ct al.l994; Katz 



(e.g. [White fc Frcnk 199l|; |CoIc 199l| ; jKatz 1992 



Navarro & 



a soluti ju Lu IhLi piubltm. LULig_y iblLajLd hum ilan iu LIil 



course of their evolution would act as negative feedback on 
the gas, limiting its cooling rate and associated star forma- 
tion. 



Primack 1999; Thacker fc Couchman 2000 Cole et al. 2000) 



Weinberg, fc Hernquist 199(j; [yepes et al. 1997 ; Somerville & 



email 



and on the chemical enrichment of the intergalactic medium 
at high redshift (e.g. [Ellison et al. 1999| ; [Schaye et al. 2000[ ). 
Efforts to detect obs ervational sig n atures of feedback are 
also being made (e.g. 



S . T . KayiQsussex . acuk 



Martin 199£ 



2001; Pettini et al. 2001) 



Theuns, Mo & Schaye 
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Recent modelling of galaxy formation, including the ef 
fects of feedback, has made use of two related techniques 
semi-analytic modelling and direct numerical simulation 
The first is a direct des cenda nt of the methodology devel- 
oped by White & Rees (1978) and provides a computation- 
ally efficient way to calculate a detailed model of galaxy for- 
mation. The merger history of dark matte r haloes is followed 
usin g ei ther sta tistic a l tec h niques (e 
1994r3ole et al. 



1994 



1997 



Kauffmann et al 
2000|: ^omerville fc Primack 1999|)"or 



A'^-b odv simulations (e.g . Kauffmann, Nusscr fc Stcinmetz 



E enson et al. 2001 ; Somerville et al. 200C| ) , the evolu- 



tion of gas is treated by means of a spherically symmetric 
cooling flow model, star formation and the effects of feedback 
are parameterised according to simple, observationally moti- 
vated rules and stellar evolutionary effects are incorporated 
by means of a stellar population synthesis model. Chem- 
ical evolution and dust extincti on or emission are readily 
incorporated in this scheme (e.g. Granato et al. 2000). Semi- 
analytic modelling has proved to be extremely successful for 
the interpretation of observational data at both low and high 
redshift, as reviewed recently by Baugh et al. (2001). 

In this paper, we are concerned with the treatment of 
star formation and feedback in gasdynamical simulations of 
galaxy formation. Numerical simulations have the advantage 
that some of the assumptions made in the semi-analytic 
treatment of gas dynamics can be relaxed. In particular, no 
artificial symmetries are imposed on a situation which is in- 
herently asymmetric. However, the application of simulation 
techniques to galaxy formation is still in its infancy because 
it places strenuous demands on computational resources, 
consequently limiting the dynamic range of resolved struc- 
tures. Early work, while making brave first attempts at simu- 
lating cosmological galaxy formation, was ha rnpcred to varv 



ing degrees by relatively poor resolution ( e.g. Ccri fc Ostrikci 



1992|; |I< atz, Hernquist fc Weinberg 1993: lEvrard, Summers 



Davi i 1994 ; Katz, Weinberg, fc Hernquist 1996 ~ |FYenk et 



al. 1996). It is only over the last few years that parallel super- 
computers have enabled simulations large enough to resolve 
present-day bright gal axies adequately within statisticallv- 



relevant volumes ( e.g. Pearce et al. 199£ ; Fardal et al. 2000 



Pearce et al. 200C). Detailed tests of varying assumptions 
and parameter values in simulations based on the Smooth 
Par ticle Hydrodynamics technique are presented in Kay et 
al. ( ^00C| , hereafter K2000) 

It is unfeasible to consider an explicit treatment of star 
formation and feedback within a cosmological simulation 
because these processes occur on much smaller scales than 
can be resolved with present techniques. The only approach 
(which is currently also used in semi-analytic models) is to 
develop phenomenological prescriptions and test their valid- 
ity by comparing their effects to as many observables as pos- 
sible. The main aim of this paper is to compare the outcome 
of adopting various star formation and feedback prescrip- 
tions in a series of small test simulations, similar to those 
analysed by K2000. We assume that feedback is produced 
only as a result of energy being injected into gas through su- 
pernovae occurring in the star forming regions that develop 
in the simulations. We also investigate the effects of varying 
the numerical resolution, albeit over a small dynamic range. 
Our primary concern is the effect of feedback on the prop- 
erties of the resolved galaxy population, mainly at z; = 0, 
but we also examine the properties of reheated gas, since 



it is these predictions that may provide some of the most 
stringent observational tests. 

We begin with a short description of our numerical 



methods in Section 
models in Section 



. We compare various star formation 
before investigating feedback in Sec- 
tion ^, focusing on both the galaxies and the distribution of 
reheated gas. The sensitivity of our results to numerical res- 
olution is explored in Section |5| Finally, we draw conclusions 
in Section ^ 



2 NUMERICAL METHOD 
2.1 The code 

The simulations were carried out using the hydra code 
^(Couchman, Thomas & Pearce 1995). This combines an 
adaptive particle-particle /particle- mesh (AP'^M ) algo rithm 
to calculate gravitational forces (Couchman 1991) with 



Smoothed Parti cle Hydrodynamics (SPH; e.g. Benz 199C 
Monaghan 1992| and references therein) to estimate hydro- 
dynamical quantities. The SPH implementation di ffers f rom 
that described in Couchman, Thoma s fc P earce (1995), as 



discussed by Thacker & Couchman (200C). In addition to 
these changes, the SPH neighbour search algorithm in hy- 
dra has been modified in order to avoid setting an artificial 
lower limit on the gas density, present in older versions of 
the code (which used a single chaining mesh to find neigh- 
bours for both the PP gravity and SPH calculations). In the 
new version, a separate neighbour list for the gas particles is 
constructed usin g the algorithm described in Appendix A2 



of Theuns et al. (199^ 



Radiative cooling is mode lled in the manner described 



by Thomas fc Couchman (1992), except that tabulated cool- 



ing rates (as a functio n of te mperature and metallicity) from 



Sutherland & Dopita (1995) are implemented. This assumes 
that the gas is in collisional ionization equilibrium and is ap- 
plicable only to radiative processes at temperatures above 
10*K . Contributions to the cooling rate at temperatures be- 
low 10*K are ignored, as are the radiative heating contri- 
butions from a photo-ionizing background. The former are 
only important in structures on scales significantly below 
the resolution of our simulations. Although photo-heating 
is important in determining the structure of the low den- 
sity gas responsible for the Ly-a forest, previous studies 
have shown that the process is unimportant for galaxies 





Quinn, Katz fc Efstathiou 199(:; 


Navarro fc Steinmetz 


19971; Weinberg, Hernquist & Katz 1997), as is the case in 



our simulations. 



2.2 Cosmological models 

We carried out simulations using two cold dark matter mod- 
els: SCDM and ACDM. The former assumes a critical den- 
sity universe, f2o = 1 and Ao — 0, Hubble constant, in units 
of lOOkms"^ Mpc~^, h = 0.5 and baryon density, Qi, = 0.06. 
The ACDM model assumes a flat, low density universe, with 



t HYDRA is publicly avail able and can be downloaded from 
tittp : //hydra. mcmaster . ca 
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Table 1. Values of key parameters of the cosmological models 
studied in this paper 



Model 


Ho 


Ao 




h 


r 


o"8 


Zi 


SCDM 


1 





0.06 


0.5 


0.5 


0.6 


24 


ACDM 


0.35 


0.65 


0.038 


0.71 


0.21 


0.9 


49 



Table 2. Summary of the star formation methods studied 
Simulation Details 



PEARCEl 

PEARCE2 

SUMMERSl 

SUMMERS2 

NAVARRO 

KATZ 

MIHOS 



5 > 5000 ; T < lO^K 
5 > 50,000 ; T < lO^K 
5 > 10 ; p > 10^25 gj.j^- 
5 > 10 ; p > 10^26 gj.^^- 
V.v < 0; p > 7 X 10-26 



V.v < 
V.v < 



■3 ; T < 3 X lO^K 
'3 ; T < 3 X lO^'K 
q; cm"^ 
-25 , 



<5 > 10; p> 10-25 gcm-3; rj > rjyn 
<5 > 10 ; Csfr = 5 X lO-"* ; Af = 3/2 



Qq — 0.35 and Aq — 0.65, Hubble constant, h — 0.71, 
and baryon density, Qt ~ 0.038. Current determinations 
of cosmological parameters (from, for example, microwave 
background anisotropics, high redshift supernovae, and the 
baryon fraction in galaxy clusters) favour the ACDM model. 
However, we do not expect our conclusions to differ signifi- 
cantly in other hierarchical models. 

Initial density fluctuations were calculated for the 
SCDM model using the appro ximate transfer function given 
by Bond & Efstathiou (1984) with shape parameter, F — 
^loh = 0.5 and amplitude set by erg — 0.6 (as required for 
agr eement with the present day abundance of rich clusters; 
e.g.|Vianna fc Liddle 199^; |Eke, Cole fc Frenk 1996]). F or th e 



ACDM model, the transfer function of Bardeen et al. (1986) 
was adopted an d F was set to F — floh exTp{—Qi,/Qo(^lo + 
V2h)) = 0.21 ( ^ugiyama 1995[ ). The amplitude of density 
fluctuations was normalised by setting erg = 0.9, as required 
by the cluster abundance argument in this case. 

The gas particle mass was held fixed for all simulations 
within each cosmological model, at ~ 5 x 10^ Mq in the 
case of SCDM and to ~ 3xlO*/i"^M0 in the case of ACDM. 
Both models assumed an initial (equivalent Plummet) soft- 
ening length of ~ 20/i~^kpc which was fixed in comoving 
coordinates until z ~ 1, after which it was fixed in physical 
coordinates to lO/i-^kpc until z = 0. The SCDM simula- 
tions were started at a redshift of 24 and the ACDM simu- 
lations at a redshift of 49 and both sets were run to z = 0. 



3 STAR FORMATION 

In this section we concentrate on comparing several star for- 
mation models drawn from the literature. We present results 
for SCDM simulations with 32768 gas and 32768 dark mat- 
ter particles, initially perturbed from a regular grid, within 
a 10 Mpc comoving box. The initial conditions and pa- 
rameters are identical to those of the fiducial simulation 
of K2000, where further details may be found. Although 
the box-size is too small to allow statistically robust results 
for the galaxy population, the simulations are only used for 
comparative purposes. 



3.1 Star formation models 

We now summarise each star formation model studied. In 
all cases, the basic procedure is the same: at the end of each 
timestep, a subset of gas particles is identified as being eligi- 
ble to form stars. A fraction of them are then converted into 
stars, which are subsequently subject only to gravitational 
interactions. Details of each method are given in Table |2| 

3.1.1 The Pearce method 



This is based on the method used by Pearce (1998). Po- 
tential star particles are identified with gas particles having 
overdensity, S > 550 and temperature, T < lO^K. We car- 
ried out two simulations, both using the same temperature 
threshold as Pearce but one with an overdensity threshold 
of 5000 (labelled pearce 1) and the other with a thresh- 
old of 50,000 (labelled PEARCE2). Our chosen overdensity 
thresholds are higher than that used by Pearce, who se- 
lected a value in order to ensure that a certain fraction of gas 
was converted into stars by ^ = 0. However, Pearce's value 
was also close to the maximum resolved SPH gas overden- 
sity in his simulations. This limit is estimated by consider- 
ing a sphere of radius 2/iniin containing A'^sph neighbours. 
(The minimum SPH smoothing length, /imin is set equal 
to 1.17 times the equivalent Plummer softening length.) 
For our choice of parameters (2/imin = 0.0234 ft-^ Mpc and 
iVsPH = 32), this limit is 5 ~ 20,000 at z = 0. We delib- 
erately set our choice of overdensity thresholds to straddle 
this resolution limit. 



3.1.2 The Summers method 



In the method used by Summers (1993) gas particles are con- 
verted into stars if they satisfy a temperature criterion, T < 
3 X lO'^K and a physical density criterion, nu > 0.1cm-''. 
Additionally, eligible gas particles must be within overdense 
regions, 5 > 10, so that star formation does not occur within 
shallow potential wells at high redshift. We performed two 
simulations using the Summers method, labelled SUMMERSl 
and SUMMERS2 respectively. For SUMMERSl, we used a den- 
sity criterion, p > lO'^^^gcm-^ (p = mHnu/X; we assume 
X = 0.76), approximately a factor of 2 lower than used 
by Summers. This threshold corresponds to a baryon over- 



density, 5 ~ 3.5 X 10^ at 



which is over an order of 



magnitude higher than our resolution limit. Because of this, 
we carried out a second simulation (summers2) for which 
we chose a density threshold of p > lO^^'' gem"''. 

3.1.3 The Navarro & White method 



We have also studied the method of Navarro & White ( 199J ) 
who applied it to follow star formation in disk galaxies. This 
is a popular method that has since been used by several 



other authors, with minor modifications (e.g. 


Steinmetz & 


Miillcr 1995; Tisscra, Lambas, & Abadi 1997 


; [Carraro et 


al. 1998; |Buonomo et al. 2000|). To summarise, regions of star 



formation are identified in which the gas is in a convergent 
flow (V.v < 0) and has a density, p > 7 x lO^^^gcm"''. 
Gas particles satisfying these criteria are split into equal 
mass gas and star particles, after a time interval. At, = 
ln2Tdyn, where rdyn oc l/^p is the local dynamical time 
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of the gas. When the gas mass is less than 5 per cent of its 
original value, the remaining fraction is converted into stars. 
We label the simulation using this method NAVARRO. 



3.1.4 The Katz, Weinberg & Hernquist method 

The fourth method investigated is that used in t he cos mo- 
logic al simulatio ns of Katz, Weinberg & Hernquist (199C ; see 
also Katz 1992). Four criteria are considered to decide if a 
gas particle is eligible to form stars: the divergence criterion, 
as used in the Navarro & White method; the physical den- 
sity threshold used as in the Summers method; an overden- 
sity threshold, 5 > 55.7, and a Jeans instability condition, 
T} > Tdyn, where rj = ft/cs is a "Jeans" timescale, h is the 
particle's SPH smoothing length and Cs = (7(7 — 
2 is the adiabatic sound-speed of the gas. (Note that this 
condition depends on the resolution of the simulation since 
h is limited by the gravitational softening.) If a gas parti- 
cle is eligible to form stars, a third of its mass is converted 
if r < 1 — exp(— At/min(rdyii, Tcooi)), where r is a random 
number between and 1, At is the timestep and TcooI oc 1/p 
is the cooling time. This method gives each particle a dual 
identity (we refer to it as a schizophrenic particle), allowing 
it to carry both a gas mass and a stellar mass. Gravitational 
forces are calculated using the total mass and hydrodynam- 
ical forces using only the gas mass. Once the gas mass has 
dropped below 5 per cent of its original value, the particle 
becomes a purely collisionless star particle. We refer to the 
simulation using this method as katz. 

3.1.5 The Mihos & Hernquist method 

The final model s tudie d is based on the method use d by Mi- 
hos & Ifernq uist ( L994 ). This uses the Schmidt law ( Schmidt 



Table 3. General properties of the simulations with star forma- 
tion at 2 = 



195£ ; see also Kennicutt 199J and references therein), an em- 
pirical (power-law) relation between the star formation rate 
per unit area and the surface density of Hi gas in late-type 
galaxies. Mihos & Hernquist convert this to a 3-dimensional 
relation, such that the fractional change in mass converted 



into stars during a timestep is AM*/Mg; 



fr P 



-'At. 



This equation is applied to gas particles with 5 > 10 and 
V.v < 0, although rather than converting part of a parti- 
cle's mass into stars, we convert all of it if AMt/Mgas > r, 
where r is a random number between and 1. We study a 
simulation using this method, labelled MIHOS, in which we 
have set A'^ — 3/2, consistent with the observed value, and 



Csb = 5 X 10" 



-1/2, 



3.2 Results 

3.2.1 General properties at z — 

Table ^ lists general properties of the 7 simulations at z — 0, 
and of the equivalent simulation without star formation (la- 
belled gas). Column 2 gives the CPU time required to com- 
plete each run, relative to GAS. (To provide a fair compar- 

i u = kT/{"/ — 1) ^mii; we take the ratio of specific heats at 
constant volume and constant pressure for a monatomic ideal 
gas to be 7 = 5/3, and the mean molecular weight of a gas of 
primordial composition to he /i = 0.59 



Simulation 


CPU 


/"cool 


/* 




GAS 


1.00 


0.35 


0.00 


74 


PEARCBl 


0.32 


0.33 


0.28 


64 


PEARCE2 


0.51 


0.34 


0.25 


70 


SUMMERSl 


0.56 


0.34 


0.16 


72 


SUMMERS2 


0.37 


0.34 


0.26 


70 


NAVARRO 


1.33 


0.35 


0.19 


69 


KATZ 


0.67 


0.35 


0.15 


71 


MIHOS 


0.48 


0.35 


0.24 


64 



ison, differences in the machine-dependent clock rate were 
removed by dividing the total time taken for each simula- 
tion by the average time per step over the first 10 steps.) 
The least efflcient method is NAVARRO, which takes around 
a third longer to complete than GAS because of the increased 
number of particles produced by star formation. The rest of 
the simulations all take around half of the time to complete 
than GAS, with the shortest being pearceI. Including a pre- 
scription for star formation reduces the CPU requirement 
because of the reduced SPH workload. 

Column 3 gives the fraction of baryons at z = in stars 
or gas with p > 500 (p) and T < 30, OOOK , labelled /cool. 
Gas that satisfies these criteria (i.e. it is cold and dense) 
is readily associated with galaxies, since the two thresholds 
enclose a region of the p~T plane containing material which 
has radiatively cooled (see, for example, §3.2 of K2000). For 
all simulat ions, including GAS, /g = 0.34 ± 0.01. As found 
by Pearce ( 1998 ) , the conversion of cold dense gas into stars 



does not significantly affect the amount of gas that cools. 

Column 4 gives the total fraction of baryons in stars, /«, 
for each run. The range of /, is much wider than for /cool, 
varying from 15 per cent (katz) to 28 per cent (pearceI). 
Increasing the overdensity threshold from 5000 ( pearceI) 
to 50,000 (PEARCE2) decreased /* by only 0.03. This smaU 
difference is due to the extra time it takes for cooled mate- 
rial to collapse to a higher overdensity before forming stars. 
SuMMERS2 produced a similar fraction of stars as pearceI 
and PEARCE2, but 40% more than SUMMERSl. The differ- 
ence in the SUMMERS simulations is also due to the choice 
of density threshold which, at z = 0, is about an order of 
magnitude higher for SUMMERSl than for pearce2. katz 
and NAVARRO produce similar values of /, (to within 20%) 
to SUMMERSl because they use the same density threshold. 
The number of stars forming at any given time in MIHOS 
is controlled by the choice of normalisation constant, Csh. 
The value, Csfr = 5 x 10"" cm^^^g'^/^s'^ was chosen (by 
trial and error) to give a similar fraction of stars at 2: = as 

PEARCE2. 

The final column in Table |^ lists the number of galaxies 
in each simulation. T hese were identifi ed using a friends-of- 
friends group-finder (Davis et al. 1985) on the cold gas and 
star particles (i.e. on those used to calculate /cooi above), 
with a linking length, 6 = 0.1 (comparable to the gravita- 
tional softening length.) Any objects with fewer star parti- 
cles than A^spH (32 for the simulations studied in this paper) 
were discarded. All runs contain ~ 70 galaxies at 2 = to 
within 10%. 
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Figure 1. Cumulative fraction of baryons in stars, normalised 
to the 2 = value, plotted against time. Also shown is the nor- 
malised fraction of cooled material (defined in the text) for the 
simulation without star formation. Rcdshift intervals are marked 
along the top of the figure. 



3.2.2 Star formation histories 

Fig. illustrates the star formation histories in the simula- 
tions by means of the cumulative fraction of baryons turned 
into stars as function of time, relative to the z = frac- 
tion (the corresponding redshift intervals are marked along 
the top of the figure) . The cumulative fraction of cooled gas 
(p > 500 (p) and T < 30, OOOK ) for the simulation with- 
out star formation, is also plotted. The instantaneous star 
formation rate is proportional to the gradient of each star 
fraction curve. 

All simulations started cooling material and forming 
stars within the first 2 Gyr [z > 2). The exact time when 
the first stars form depends upon two factors: the resolu- 
tion of the simulation, which determines the formation time 
of the first resolved dark matter haloes, and the range of 
gas densities where star formation occurs. Since the resolu- 
tion is fixed in all our simulations, only the second factor is 
relevant. 

The star formation history in pearceI is similar to that 
of cooled material, i.e. it rises sharply from z = 5 until z ~ 2, 
before rolling over, implying a decreasing star formation rate 
at lower redshift. However, in pearce2, the delay caused 
by a higher overdensity threshold causes star formation to 
begin later and to proceed at a slower rate until around 5 
Gyr {z ~ 1), after which the z = normalised star formation 
rate is higher than in pearceI. 

SUMMERS 1, NAVARRO and KATZ have the same histories 
at all times, demonstrating that the additional features of 
the latter two models are irrelevant in determining the aver- 
age star formation rate within the simulations. (We exclude 
the latter two simulations from further discussion since they 
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produce nearly identical results to summers for all quanti- 
ties studied in this paper.) 

These curves are very similar to those from the pearce 
simulations until z — 1 — 2, where a "kink" develops, due to 
the behaviour of the gravitational softening. At early times, 
the softening is constant in comoving coordinates, and so 
the minimum resolved length scale in physical coordinates 
grows in proportion to the expansion factor, a = + z), 
leading to a decrease (as (1 + z)^) in the maximum resolved 
physical density. By z ~ 1.7, this density has dropped to 
p ~ 10~^®gcm~^, the summers 1 density threshold. This 
suppresses further star formation until z = 1, when the soft- 
ening becomes constant in physical coordinates, allowing the 
star formation rate to increase again. SUMMERS2 exhibits a 
similar star formation history to pearceI. In these cases no 
kink is evident since the density thresholds are always lower 
than the maximum resolved density. 

Finally, in mihos, the relative fraction of stars rises more 
sharply than other curves at high redshift but fiattens off 
by 2 = 0. For this run, an approximate minimum density 
threshold can be calculated by first equating the relevant 
star forma tion timescale, Tsb = MgasAt/AA'/« — l/CsfrPgas 
(c.f. ^3.1.5| ) to the lookback time for an f2 = 1 universe, 
Tiook = (2/3ii/'o)(l — (1 + ^)~^^^)', at higher densities, Tsh < 
Tiook- Assuming that 5 ^ p/ (p) (i.e. for S » 1), the min- 
imum overdensity scales as Smin oc ((1 -I- z)^^'^ — 1)~^, a 
decreasing function of redshift. At high redshift, this over- 
density threshold was sufficiently low to allow efficient star 
formation to take place. However, at lower redshift, the de- 
crease in the density of the Universe and in the available 
time left to form stars before 2 = 0, causes the number of 
gas particles above this threshold to drop rapidly. 

3.2.3 Cold gas fractions 

In Fig. ^ we plot the median mass ratio of cold gas to stars 
against the stellar mass of each galaxy with 32 or more 
star particles and a non-zero gas mass, in bins of width, 
AlogAf* = 0.35. Error bars represent 10 and 90 percentile 
points of the distribution in each bin. There is a clear trend 
for more massive galaxies to contain proportionally less cold 
gas. This trend refiects the distribution of mean stellar ages 
with galaxy mass: the more massive galaxies in all star for- 
mation models harbour, on average, older stellar populations 
than less massive galaxies. Therefore, cooled baryons have 
had more time to collapse and cross the imposed threshold 
for star formation. 

The cold gas to star mass fraction in the brightest galax- 
ies (above M« = W^^h~^ M©, approximately the mass of an 
L* galaxy) is less sensitive to the choice of star formation 
prescription than for less massive galaxies. However, below 
M» = 10^^h~^ Mq, simulations with a higher density thresh- 
old (e.g. SUMMERS2) contain significantly more cold gas than 
simulations with a lower threshold. For example, the median 
ratio of cold gas to stars at M* = 2.4 x 10^°h~^ Mq is ~ 0.4 
for SUMMERS2 but is 1.6 for SUMMERSI. Again, this trend 
refiects the reduction in star formation activity when the 
density threshold is high. 

Observations of atomic and molecular gas in disk galax- 
ies can, in principle, be used to place constraints on plausi- 
ble star formation models. However, this can only be done 
in a rather crude manner because the simulations do not 
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Figure 2. The median ratio of cold gas to stars at 2 = 0, plotted 
against tlie stellar mass of each galaxy with non-zero gas content, 
for simulations with star formation. Error bars represent the 10 
and 90 percentile points of the distribution in each bin. 



Figure 3. Cumulative mass functions for galaxies in simulations 
with star formation at z = 0. The galaxy mass includes cold 
gas and stars. For comparison, the mass function for GAS is also 
plotted. 



have enough resolution to distinguish between morphologi- 
cal types. By selecting those galaxies with cold gas we are 
at least able to identify candidate disk galaxies in the simu- 
lations. From the data plotted in Fig. 9 of Cole et al. (2000), 
L* disk galaxies [M/Lb ~ 2 and Mg ~ —19.5) have 
Mcg/M* ~ 0.1 — 0.3, values which do not increase signif- 
icantly for galaxies an order of magnitude less massive. This 
is broadly consistent with our results for simulations with 
low density thresholds, but inconsistent with those with high 
thresholds (extrapolating the result for SUMMERSl down to 
M« = lQ^°h~^ Mq implies a ratio of at least 2). 



3.2.4 Galaxy mass and luminosity functions 

Cumulative mass functions for the galaxy population at 
z — are plotted in Fig. ^ For comparison, the mass func- 
tion for simulation GAS is also displayed. All the mass func- 
tions are very similar, demonstrating that the abundance of 
galaxies over the mass range resolved in these simulations is 
insensitive to the inclusion of star formation. 

We can assign luminosities to the galaxies in the simu- 
lations using a stellar population synthesis model. We cal- 
culate ii'-band luminosities since the hght in this band is 
less affected by dust than at other bands and, since it is 
dominated by late-type stars, it provides a more robust es- 
timate of the underlying stellar mass of the galaxy. To cal- 
culate galaxy luminosities, the formation time, i/, of each 
star particle, m,, was used to calculate the luminosity per 
unit mass at the present day, l{to — tf). The total magnitude 
is then 
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Figure 4. Cumulative ii"— band luminosity functions of the galax- 
ies at 2 = for simulations with star formation. The Schechter 
fit to the luminosity function from the 2dF and 2-Mass surveys 
(Cole et al. 2001) is also plotted (solid curve), with the value of 
shown as a vertical arrow. The horizontal arrow indicates 
the shift that would be required for simulation PEARCEl to agree 
with the observational data at 



Mk = 51og/i + Mk{0) - 2.5 log ^ l{to 



■tf)m,, 



(1) 
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where M/f (0) is the zero-point and A^, is the number of star 
particles in the galaxy. The function, l{tQ~tf) was generated 
using the stella r population synthesis models of Bruzual & 
Chariot ( [L993| ), assuming so lar m etallicity and the stella r 
IMF proposed by Kennicutt ( |l983| see also |Cole et al. 2000| ). 
The total magnitude of each galaxy can be rescaled to take 
into account non-visible stellar material (i.e. a population 
of brown dwarfs) by adding the term 2.51ogT, where T = 
1 + M(brown dwarfs) /M (visible stars) (|Cole et al. 2000[) . 



Cumulative luminosity functions for galaxies in our sim- 
ulations at z = are shown in Fig. |^. Only galaxies with 
32 or more star particles (or equivalent) are included. The 
best-fit Schechter function to the _K'-band galaxy luminosity 
function from the 2dF and 2- Mass surveys is also plotted as 
a solid curve (Cole et al. 2001), with the value of marked 
with a vertical arrow. The model luminosity functions are 
plotted assuming T = 1; the shift required for simulation 
pearceI to agree with the data at M'^ is indicated by a 
horizontal arrow. This implies a value of T = 3, i.e. twice 
as much mass is required in brown dwarfs than in visible 
stars, which is in c lear contradiction with observations (see 
Balogh et al. 200l[). 



3.3 Summary 

We havp examined \farious prescriptions in the literature for 
qs into stars in roHmo1np;ica1 simirlatinns The de- 



turning 



terminant factor in each model is the range of density (usu- 
ally defined by a lower limit, or threshold) at which star 
formation is allowed to proceed efficiently: simulations with 
lower thresholds naturally form more stars than simulations 
with higher thresholds. The cooling rate of the gas is largely 
unaffected by the inclusion of a star formation mechanism. 
In these models, the star formation efficiency determines 
the fraction of galactic baryons that end up in the form of 
cold gas or stars. Simulations with a high density thresh- 
old contain more cold gas in (particularly sub-L*) galaxies 
than simulations with low density thresholds; the former are 
discrepant with observational data. All of the models over- 
predict the abundance of galaxies at a given magnitude, or 
equivalently produce galaxies that are (about 0.5-1 magni- 
tudes) too bright at a given abundance. In the next section, 
we consider feedback mechanisms that act to reduce galaxy 
luminosities. 



4 FEEDBACK 

In this section we study the effects of energy injection asso- 
ciated with Type II supernovae explosions in a cosmological 
simulation. On energetic grounds alone, these events must 
be important during galaxy formation, and their role as a 
negative feedback mechanism allows ab-imtio hierarchical 
models to match th e observed galaxy luminosity fun ction 
(IWhite & Rees 197^; IWhite & Frenk 199l|; ICole 199l|). Ide- 



ally, we would like to be able to resolve individual supernova 
events and model the multi-phase intergalactic medium in 
our simulations. However, this is impractical, both because 
of computational limitations and because of the lack of a de- 
tailed understanding of the processes involved. At present, 
the only practical strategy is to include the effects of super- 
nova feedback phenomenologically. 



Phenomenological feedback models within numerical 
simulations have already been considered by a number of 
authors studying galaxy formation/evolution. For example, 
supernova fee dback was first included within SPH simula- 
tions by Katz(1992), who studied the collapse of a galaxy- 
sized density perturbation, and used a scheme where every 
star particle supplied neighbouring gas particles with ther- 
mal energy (heat). However, this had very little effect on the 
galaxy's properties, because the high gas density at the loca- 
tions where energy is injected implies a cooling time which 
is much shorter than the dynamical time. As a result, the 
energy injected is simply radiated away. For future reference 
we will refer to this scheme for energy injection as "thermal 
feedback" . 

An alternativ e fee dback scheme was proposed by 
Navarro & White (1993). In this model, a fraction, /„, of 
supernovae energy is supplied to neighbouring gas particles 
in the form of kinetic energy, with the rest being added as 
heat. Like Katz, Navarro & White found that for low values 
of /„ (i.e. almost all thermal energy), very little feedback oc- 
curs and consequently the star formation rate is very high. 
However, for higher values of /„ the star formation rate is 
reduced because kinetic perturbations are able partially to 
halt the local collapse of a density perturbation. We refer 
to this method as "kinetic feedback". Kinetic feedback has 



since been used by other authors (e.g. Mihos fc Hernquist 



1994) 



Subsequent attempts have been made to boost the ef- 
fects of thermal feedback by preventing the gas from im- 
mediately re-radiating th e supe rnova energy. This was first 
considered by Gerritsen (1997) who argued that previous 
feedback schemes did not produce an overall disk morphol- 
ogy resembling those observed, with spherical cavities, or 
bubbles, filled with hot gas. They found that injecting all the 
energy into one particle as heat and preventing it from cool- 
ing for a finite length of time, had the de sired effect. A more 
elegant approach was taken by Springel (2000), who created 



a second pressure term in the equations of motion, represent- 
ing turbulent pressure caused by supernovae. Yet an other 
example is the study by Thacker & Couchman (200C), who 



considered a scheme in which a "cooling" density for each 
gas particle injected with energy is obtained from a simple 
pressure equilibrium argument. This density is lower than 
the standard SPH density, and so lengthens the particle's 
cooling time, before resorting back to the SPH density after 
a finite period. 

Attempts to simulate galaxy formation in a large region, 
starting from realistic cosmological initial conditions, and 
including a model for feedback have been made by severa l 



groups using Eulerian gr id codes (e.g. Cen fc Ostriker 1992 



199S|; [Yepes et al. 1997[) and SPH codes ( Rummers 199; 



Katz, Weinberg, & Hernquist 1996). The former technique 



is hindered by poor spatial resolution, with galaxies being 
several factors smaller than the minimum resolved scales (al- 
though this situation is improving, particularly with adap- 
tive mesh codes). Studies using SPH simulations have fo- 
cussed on thermal feedback, and like Katz (1992), have also 



found that adding thermal energy to the gas makes little 
difference to the resulting galaxy population. In the next 
subsection we will describe our own implementations of ther- 
mal and kinetic feedback in our simulations, which we r egard 
as a natural follow-up to the work of Summers (1993) and 
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Katz et al. ( |l996| ). As we shall see, our models are capable of 
producing strong feedback, reducing the amount of material 
which cools and forms stars, and resulting in a more realistic 
A'-band luminosity function. 



4.1 Simulation details 

4-1.1 Cooling and Star formation 

We have studied the effects of feedback in the ACDM 
and SCDM cosmological models described in § 2.2. Simi- 



lar trends with respect to variation of the parameters are 
present in both models, so we will focus mainly on ACDM. 
We impleme nt th e "decoupling" technique described by 
Pearce et al. (200C) whereby gas particles with temperature 
below T — 12, OOOK are ignored when estimating the gas 
density of hot particles with temperature above T = 10"" K . 
This technique effectively prevents the unphysical cooling of 
hot gas induced purely by the presence of nearby cold, dense 
galactic gas. (Note, however, that this problem is largely 
circumvented due to the inclusion of a star formation pre- 
scription which removes the majority of galactic material 
from the SPH calculations.) This modification is a crude at- 
tempt to model a multiphase gas with SPIf and has b een 
disc ussed further elsewhere (e .g. K2000; Croft et al. 200C ; see 
also [Ritchie fc Thomas 200 1| ). 

Star for mation was incorporated using the Pearce 
method (c.f. § 3.1.1), assuming an overdensity threshold of 



5000 and a temperature threshold of 30, OOOK (essentially 
the pearceI criteria). When a gas particle satisfies both 
criteria it is instantaneously converted into a star particle. 
In the ACDM simulations we assumed an unevolving gas 
metallicity of Z = OSZq while for the SCDM simulations, 
we retained the metallicity of Z = 0-5Zq assumed in the pre- 
vious section. Note that our choice of metallicity results in 
an overestimate of the cooling rate of gas at high redshift, 
where Z is considerably smaller. Self-consistent models of 
feedback and metal enrichment within simulations will be 
presented in a future paper. 



4-1-2 Supernova energetics 

We adopt the same parameter values for the supernovae 
energetics as Katz (1992) and Katz et al. ( 1996| ). To sum- 
marise, the amount of energy available from supernovae is 
calculated assuming that stars with mass above SM© release 



lO^rg of energy into the surrounding gas. A Miller & Scalo 



(1979) IMF is assumed, with lower and upper mass limits of 
O.IM© and IOOMq respectively, implying a specific energy 
released in supernovae of ~ 3.7 x 10^^ ergg~^, equivalent to 
a temperature of 1.8 x lO'^K . The size of the timestep in the 
simulations (typically 3Myr) is about one order of magni- 
tude smaller than the lifetime of an 8M0 star (t8 ~ 20Myr). 
For this reason, the energy was released gradually, using the 
following function 



usn{t) = 3.7 X 10'^ exp {-{t - to)/Ts) Tg ' erg 



(2) 



where to is the time when the star formation event occurred. 
Energy was released until t ~ to = 200Myr, when -usn be- 
comes negligible. 



4-1-3 Feedback models 

Our prescription for thermal feedback works as follows. 
When a gas particle is to be converted into a star particle, 
it first passes through an intermediate phase, in which it is 
labelled a supernova particle. Supernova particles feel the 
gravitational force only, but the SPH algorithm maintains 
both their densities and neighbour list. Smoothing lengths 
are fixed at the minimum value, ftmin, and each gas neigh- 
bour within the SPH smoothing radius (2/imin) receives a 
proportion of the thermal energy available at that time. For 
supernova particle i, the energy given to gas particle j, over 
the timestep At, is 



H{2hmin - |ri - Fjl) 



J2"2T H{2h^,^-\r,-r,\) 



(3) 



where < e < 1 is the feedback efficiency parameter, and 
H is the Heaviside step function. We choose not to weight 
the contributions using the SPH kernel since this as an un- 
necessary complication in cosmological simulations, where 
the forces within galaxies are softened. Furthermore, the ef- 
ficiency of energy transfer will vary depending on the num- 
ber of neighbours of each supernova particle. Our method 
ensures that the same amount of energy per supernova is al- 
ways transferred to the gas and so the efficiency is controlled 
solely by the choice of e. 

As discussed above, this scheme is expected to have 
little effect on the galaxy properties, since the gas will re- 
radiate the energy in a short time. To cir cumv ent this prob- 
lem, we follow the method of Gerritsen (1997) and prevent 
the gas from cooling over a period of time. This timescale 
is controlled using the parameter r, which is the length of 
the non-radiative period in units of the age of the universe 
at z = 0. 

We also study kinetic feedback, where instead of heat- 
ing the gas, it receives a velocity perturbation radially away 
from the supernova particles. The total magnitude of the 
velocity added to gas particle j, from supernova particle i, 



AwsNj = , / 2 e usN.i At ■ — — 



(4) 



This scheme is the "momentum" implementation of Navarro 
& White 19931 



4.2 Results 

We present results from 6 different parameter choices for 
the thermal feedback model and 3 for the kinetic model. 
Specifically, for simulations with thermal feedback, we con- 
sider r = 0.01, 0.1, 1.0 for e — 0.1, 1.0. For simulations with 
kinetic feedback, we consider e — 0.001,0.01,0.1. Details of 
the feedback models for ACDM and SCDM simulations are 
summarised in Table ^ 



4-2-1 General properties at z — 

We identify galaxies by applying a friends-of- fr iends group- 
finder to the gas and star particles as in Section 3.2.1 , except 
that we use a linking length b = 0.2. This gives a slightly 
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Table 4. General properties at ^ = of the simulations with 
feedback 

Cosmology Feedback e t /g^i ^gal 

ACDM None N/A N/A 0.35 34 

SCDM None N/A N/A 0.26 53 

ACDM Thermal 0.1 0.01 0.28 30 

ACDM Thermal 0.1 0.1 0.17 26 

ACDM Thermal 0.1 1.0 0.096 17 

ACDM Thermal 1.0 0.01 0.11 13 

ACDM Thermal 1.0 0.1 0.083 12 

ACDM Thermal 1.0 1.0 0.065 13 

SCDM Thermal 0.1 0.01 0.19 49 

SCDM Thermal 0.1 0.1 0.099 35 

SCDM Thermal 0.1 1.0 0.061 21 

SCDM Thermal 1.0 0.01 0.076 24 

SCDM Thermal 1.0 0.1 0.047 14 

SCDM Thermal 1.0 1.0 0.037 11 

ACDM Kinetic 0.001 N/A 0.34 29 
ACDM Kinetic 0.01 N/A 0.22 20 
ACDM Kinetic 0.1 N/A 0.056 11 

SCDM Kinetic 0.001 N/A 0.23 52 
SCDM Kinetic 0.01 N/A 0.15 35 
SCDM Kinetic 0.1 N/A 0.046 11 

larger number of galaxies than b — 0.1. (In simulations with- 
out feedback, both values of b pick out the same galaxies, but 
when feedback is included a larger linking length is required 
because the galaxy material becomes slightly distended by 
the disturbance produced by the injection of supernova en- 
ergy. Since the internal structure of galaxies is unresolved, 
the exact choice of linking length is unimportant provided 
it captures the majority of galaxy material.) We also de- 
mand that each galaxy should contain a minimum of 32 
star particles, regardless of the number of gas particles, in 
order to provide a reasonable estimate of its luminosity. The 
minimum baryonic mass of a galaxy in these simulations is 
l.Ox 10^°/i"^ Mq for the ACDM model and 1.6x 10"/i"^ M© 
for the SCDM model. 

Column 5 of Table ^ gives the fraction of baryons (cold 
gas and stars) in galaxies, labelled /gai. For simulations with 
the same feedback parameters, the fraction of galaxy mate- 
rial is always higher in the ACDM case than in the SCDM 
case though the trends with varying e and r are the same. A 
larger feedback efficiency, e, causes a decrease in the fraction 
of galactic material at 2 = (for both thermal and kinetic 
feedback). In the thermal case, increasing the efficiency in- 
creases the pressure of reheated gas particles, with the result 
that more material can rise out of galaxies, before it can re- 
cool and form stars. In the kinetic case, a larger e causes 
an increase in the fraction of gas particles with velocities 
greater than the escape velocity of each galaxy. For a given 
value of e in the simulations with thermal feedback, a larger 
r results in a smaller /gai. This is because the longer the re- 
heated gas is prevented from cooling, the easier it is for this 
gas to escape from the galactic potential. (Note that for the 
largest value of t considered, t = 1, reheated gas can never 
re-radiate its energy before the end of the simulation.) 

The number of resolved galaxies in each simulation at 
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Figure 5. Cumulative fraction of baryons in stars plotted against 
time, in Gyr, for ACDM simulations without feedback (solid lines) 
and with feedback. For thermal feedback, dotted lines correspond 
to T = 0.01, short-dashed to r = 0.1 and long-dashed to t = 1; for 
kinetic feedback, the same sequence corresponds to e = 0.001, 0.01 
& 0.1. Model parameters are given in the legend in each panel. 

2; = is listed in Column 6. This number is sensitive to the 
feedback mechanism. In general, SCDM simulations contain 
more galaxies than the equivalent ACDM simulations, al- 
though the exact values depend on both the mass resolution 
of the simulation and the shape of the galaxy mass function 
near this threshold. 

^.2.2 Star formation rates 

In Fig. we illustrate the effect of feedback on the global 
star formation rate, by plotting the cumulative fraction of 
baryons in stars as a function of time for ACDM models. 
In each panel we compare results between simulations with 
the same feedback prescription, but with different values of 
one parameter. For thermal feedback, dotted lines represent 
results from simulations with t — 0.01, short-dashed with 
r — 0.1 and long-dashed with r = 1; for kinetic feedback, the 
same sequence represents results for e = 0.001,0.01 & 0.1. 
The equivalent simulation without feedback is also plotted 
as a solid line. 

Consider first the case of thermal feedback. For e — 
0.1, the efficiency of star formation is sensitive to the value 
of T (when varied over 2 orders of magnitude). However, 
increasing e by a factor of 10 produces only small differences 
in the star formation rate when r is varied. In the latter case, 
enough energy is deposited into the gas that even on the 
shortest timescale considered, r = 0.01 (of order 100 Myr), 
much of the reheated material is unable to cool and form 
stars before the end of the simulation. Kinetic feedback can 
be very effective at reducing the amount of star formation, 
even though gas above lO^K is able to cool radiatively at 
all times. Although supplying 0.1 per cent of the available 



ACDM None 

SCDM None 

ACDM Thermal 

ACDM Thermal 

ACDM Thermal 

ACDM Thermal 

ACDM Thermal 

ACDM Thermal 

SCDM Thermal 

SCDM Thermal 

SCDM Thermal 

SCDM Thermal 

SCDM Thermal 

SCDM Thermal 

ACDM Kinetic 

ACDM Kinetic 

ACDM Kinetic 

SCDM Kinetic 

SCDM Kinetic 

SCDM Kinetic 



In / A 


IVT / A 

In / A 





35 


IVT / A 

In / A 


IVT / A 

In / A 





26 


0.1 


0.01 





28 


0.1 


0.1 





17 


0.1 


1.0 





096 


i.U 


U.Ul 


U 


1 1 
ii 


1 n 
i.U 


n 1 
U. i 


U 




1.0 


1.0 





065 


n 1 
U. i 


U.Ui 


U 


1 


n 1 
U. i 


n 1 
U. i 


U 


non 

uyy 




1 n 

I.U 


n 
u 




1.0 


0.01 





076 


1.0 


0.1 





047 


1.0 


1.0 





037 


0.001 


N/A 





34 


0.01 


N/A 





22 


0.1 


N/A 





056 


0.001 


N/A 





23 


0.01 


N/A 





15 


0.1 


N/A 





046 



© 0000 RAS, MNRAS 000, 000-000 



10 



Scott T. Kay et al. 




Figure 6. Thermal history (density, temperature and entropy 
as a function of time) of a gas particle in the ACDM simulation 
with thermal feedback (e = 1, r = 1) that would otherwise have 
cooled and formed stars. The density threshold for star formation 
is plotted as a dotted line. Redshift intervals are labelled along 
the top of the figure. 



energy makes no significant difference to tlie star formation 
rate, increasing tlie efficiency to 1 and 10 per cent lowers tiie 
fraction of stars formed by 2 = to around 2/3 and 1/4 of 
the no-feedbacic values respectively. 



4-2.3 The fate of reheated gas 

Virtually all the gas that is reheated, either thermally or 
kinetically, has a temperature in the range 10* — lO^K . The 
temperature to which it is reheated and its subsequent ther- 
mal evolution depends on the specific feedback mechanism. 
For thermal feedback, at least 80% of the remaining gas that 
was once reheated is hotter than lO^K since any gas that 
can cool again is rapidly turned into stars. In the case of 
kinetic feedback, the final temperature of the reheated gas 
depends on the value of e: a greater efficiency leads to hotter 
gas. Unlike in the thermal models, gas in the kinetic mod- 
els is not reheated instantaneously but has to thermalise its 
energy, a process which will be parti cularly susceptible to 
numerical resolution (Katz et al. 1996). 



To investigate the fate of reheated gas, we first illus- 
trate the typical history of a reheated gas particle that 
would otherwise have cooled within a galaxy and been con- 
verted into a star particle. We have selected a particle in 
the ACDM simulation with thermal feedback, for e = 1 and 
r = 1. The density (n/ cm~^), specific "entropy" (defined as 
s = kT/n'^''^) and temperature of this particle are plotted 
in Fig. y. Also shown (dotted line) is the density at each 
redshift equivalent to the overdensity threshold for star for- 
mation, i.e.n > 5000 (p) //^mn. 

Initially {z > 5), the particle cools adiabatically 
(i.e. isentropically) from its initial temperature of lOOK due 



to the expansion of the Universe. (Note that our choice of 
initial temperature is arbitrary but is too low to have any 
effect on the subsequent evolution of the particle. The true 
temperature of the gas at high redshift is determined by the 
intensity of the UV background, which we ignore for this 
study.) At z = 5 the particle has fallen into a halo where it 
is shock-heated to a temperature of lO'^K : the density is suf- 
ficiently high that the particle is able to radiate any excess 
thermal energy within a timestep (K2000). At this point 
the density does not increase (as expected during shock- 
heating) since the region is not fully resolved by the SPH 
algorithm. A few of the particle's neighbours are not part 
of the overall collapse (leading to an underestimate of the 
gas density) even though the net fiow is converging, which 
results in viscous heating. Shortly after this time, however, 
the density of the particle starts increasing as it condenses to 
form part of a "galaxy" within a forming dark matter halo. 
(The halo, with Vc ~ 60kms~^, contains approximately 25 
particles and so is also poorly resolved at this stage.) The 
density reaches a peak value of 0.1 cm~"^ at z; ~ 3. However, 
before the particle crosses the threshold for star formation, 
it absorbs supernovae energy which raises its temperature 
by about 2 orders of magnitude, to lO^K . Correspondingly, 
the density decreases by 3 orders of magnitude (aided by 
the fact that since the particle is at T > lO'^K , its SPH 
density is determined only by other hot particles) and its 
entropy jumps to ~ lOOkeVcm^. Apart from a few minor 
heating events, the particle evolves approximately isentrop- 
ically for the remaining time {z < 3) as it moves into less 
dense environments. 

In order to quantify the fate of gas reheated by super- 
novae energy, we plot in Fig. ^ the fraction of reheated gas 
particles in the ACDM simulations as a function of time, 
which lie in regions of different overdensity. The solid lines 
show the entire reheated fraction, the long-dashed lines the 
fraction with 5 < 500, the short-dashed lines the fraction 
with 5 < 50 and the dotted lines the fraction with 5 < 10. 
The first cut approximately selects gas that has been ex- 
pelled from galaxies (blow-out), the intermediate overden- 
sity preferentially picks out gas expelled from haloes alto- 
gether (blow-away) , while the final criterion indicates par- 
ticles that have been ejected into weakly overdense regions 
such as those responsible for the Ly-a forest observed in 
quasar spectra. 

The first reheated particles appear at around z = 5, 
when the first stars form and the associated supernovae in- 
ject energy into the gas. The reheated fraction increases 
steadily to 2: = in all simulations except in the case of 
thermal feedback with e = 0.1 and t = 0.01. In this simu- 
lation, the amount of reheated gas that can cool again and 
form stars at z < 1 outweighs the total amount of gas that is 
reheated. As expected, at low redshift, the feedback models 
with the highest star formation rates have the lowest fraction 
of reheated gas. However, the amount of reheated material is 
lower in the kinetic feedback models at any given time, than 
in a thermal model which produces approximately the same 
fraction of stars. For example, both the thermal feedback 
model with e = 1,t = 1 and the kinetic feedback model 
with e = 0.1 predict that 5-6% of the baryons are locked 
into stars in resolved galaxies (M > W^^Mq) at 2 = 0. 
However, ~ 34% of the baryons are reheated in the thermal 
case, whereas only ~ 20% are reheated in the kinetic case. 
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Figure 7. The fraction of gas reiieated by supernovae energy as a function of time, in Gyr, lying in regions of different overdensity. Solid 
lines illustrate the total reheated fraction, long-dashed lines the fraction with 5 < 500, short-dashed lines the fraction with 5 < 50 and 
dotted lines the fraction with 5 < 10. Redshift intervals are labelled along the top of the figure. 



In all our simulations, a significant fraction of the gas 
that absorbs supernova energy escapes from galaxies and 
haloes. Kinetic feedback is more efficient at transporting the 
gas into the intergalactic medium than thermal feedback. 
In the latter models, as much as 10% of the baryons are 
reheated and expelled into weakly overdense regions {5 < 
10) by z = 0, so that the ratio of diffuse reprocessed gas 
to stars is as high as 2. Kinetic feedback models predict a 
similar ratio, even though the total fraction of reheated gas 
is much lower. 



estimate it, we first identify dark matter haloes at the simu- 
lation output times {z = 5, 3, 2, 1, 0.5, 0). Haloes are defined 
as overdense spheres enclosing an average density Ac per (z), 
where pcr{z) = 3H' ^(z) /StvG; we set Ac = 1 78 0.°-'^^ (z) for 
the ACDM model ( |Eke, Cole fc Frenk 199^ ). We then tag 



each particle with the latest output time when it was inside 
a halo, and note the corresponding halo mass. The circular 
velocity of a halo is related to its mass and redshift by 



Vc = 740 M- 



1/3 



H{z) 
Ho 



(5) 



4-2.^ Galactic winds 

One of the effects of energy injection is to expel gas from 
haloes in the form of a galactic wind. The time at which 
the gas is expelled is only approximately known in the sim- 
ulations because of the relatively sparse output times. To 



where M14 is the total mass of the halo in lQ^'^h~^ Mq. A re- 
heated gas particle is deemed to have been expelled from the 
halo when it first moves beyond the virial radius, providing 
it was identified with a resolved halo (32 or more particles) 
before the reheating event and also prior to expulsion. All 
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Figure 8. The median mass fraction of gas expelled from haloes 
as a function of the circular velocity of the halo from which it was 
expelled. Error bars illustrate 10 and 90 percentile points within 
each bin. Stars, triangles and squares represent thermal feedback 
simulations with t = 0.01,0.1 & 1.0, and kinetic feedback simu- 
lations with e = 0.001,0.01 & 0.1 respectively. 



Figure 9. The 2 = median mass ratio of cold gas to stars 
plotted against the stellar mass of each galaxy. Error bars repre- 
sent 10 and 90 percentile points within each bin. Stars, triangles 
and squares illustrate results for thermal feedback simulations 
with r = 0.01,0.1 & 1.0, and kinetic feedback simulations with 
e = 0.001,0.01 & 0.1 respectively 



expelled particles are then grouped together to calculate the 
expulsion mass for each halo. 

In Fig. 1^, we plot the median mass fraction of expelled 
gas as a function of the circular velocity of the halo from 
which the gas was blown away. (We define the halo gas 
mass as Mhaio f^b/f^o , where Mhaio is the total halo mass.) 
Error bars illustrate 10 and 90 percentile points of the dis- 
tribution within each bin. Stars, triangles and squares repre- 
sent results from ACDM thermal feedback simulations with 
T = 0.01,0.1 & 1.0, and kinetic feedback simulations with 
e = 0.001,0.01 & 0.1 respectively. Although the range of 
resolved circular velocities is limited in these simulations. 
Fig. 1^ shows a clear trend of decreasing expulsion fraction 
with increasing circular velocity for nearly all cases, with 
only a few per cent of gas escaping from Milky Way sized 
haloes (14 ~ 220kms~^) and larger. Increasing the avail- 
able energy (large e) causes a larger fraction of the gas to be 
expelled, particularly in the haloes with the lowest values 
of Vc. A larger value of e raises the mean temperature of 
the reheated gas, allowing a greater fraction to escape from 
the smaller systems. However, the average cooling times for 
these systems are considerably shorter than their dynami- 
cal times, and so the amount of gas escaping in the thermal 
feedback simulations is lower for smaller values of the delay 
timescale, r. 



4-2.5 Galaxy properties 

We now investigate the effects of feedback on the present 
day galaxy population. Firstly, we examine its effect on the 
amount of cold gas in galaxies, by plotting the mean ratio 
of cold gas mass to stellar mass, Mcg/M, , against A/, for 



galaxies (with non-zero gas content) in the ACDM simula- 
tions with 32 or more particles. The points in Fig. ^ represent 
the median values within bins of width, AlogM* = 0.5, and 
error bars represent 10 and 90 percentile points within each 
bin. The trend of decreasing cold gas fraction with galaxy 
mass is similar to that seen in Fig. 2 for simulations without 
feedback, but all our feedback models predict median cold 
gas fractions below 0.5 for galaxies with M* > 10^° Mq. 
This agrees well with the observational data summarised by 
Cole et al. ( ^000| ; see Section ^.2.3[ ) 

The overall reduction in star formation activity caused 
by feedback has a substantial impact on the galaxy luminos- 
ity function. Cumulative if -band luminosity functions are 
shown in Fig. ^ for the ACDM simulations. Luminosities 
for each galaxy (consisting of at least 32 star part icles) were 
calculated using the method described in Section 3.2.4, The 



best cumulative Schechter function fit to the data from Cole 
et al. (2001) from the 2dF and 2-Mass surveys is plotted as 
a solid curve, with the value of indicated with a vertical 
arrow. 

Results from the thermal feedback simulation with 
e — 0.1, T = 0.01 and the kinetic feedback simulation with 
e = 0.001 are similar to that from a simulation with no feed- 
back. These all produce a luminosity function which is ~ 2 
magnitudes brighter than observed. As the strength of the 
feedback is increased (by increasing e and/or r), the stellar 
mass and hence the luminosity of the galaxies are reduced. 
Both the thermal model with e = 1 & r = 1 and the kinetic 
model with e = 0.1 predict about the correct number of 
L*K galaxies, although the latter model underestimates the 
abundance of fainter galaxies. 



© 0000 RAS, MNRAS 000, 000-000 



Star formation and supernova feedback in cosmological simulations 13 



Thermal feedback 
e = 0.1, T = 0.01,0.1,1.0 



Thermal feedback 
e=1.0, T = 0.01,0.1,1.0 





Table 5. Properties of the simulations with feedback and varying 
resolution 



Kinetic feedback 
6 = 0.001,0.01,0.1 
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Figure 10. i^T-band luminosity function of galaxies in ACDM 
simulations with feedback. Dotted, short-dashed and long-dashed 
lines represent thermal feedback results with r = 0.01,0.1 Sz 1.0, 
and kinetic feedback results with e = 0.001,0.01 & 0.1 respec- 
tively. Also shown (solid curve) is the best— fit Schechter function 
to the result from the 2dF and 2-Mass surveys (Cole et al. 2001). 



4.3 Summary 

We have investigated two ways of including feedback from 
supernovae energy in cosmological simulations, in order to 
reduce the amount of gas that c ools and fo rms stars. The 
first method, thermal feedback (Katz 199i), adds thermal 
energy to surrounding gas (with efficiency controlled by the 
parameter e) and prevents such material from recooling for a 
period (controlle d by the parameter r). T he second method, 
kinetic feedback (Navarro & White 1993), adds velocity per- 



turbations to the gas radially away from the stars, and al- 
lows the gas to cool at all times. Increasing e and/or r in- 
creases the strength of feedback, reducing the global star 
formation rate, but without significantly affecting the aver- 
age ratio of cold gas to stars in galaxies. The feedback works 
in both models by reheating gas which has already cooled, 
with stronger feedback resulting in more gas escaping into 
regions of low overdensity. The efficiency with which gas is 
expelled from dark matter haloes is higher in haloes with 
Vc ~ lOOkms"^ than in haloes with higher circular veloc- 
ity. By limiting the star formation rate, feedback lowers the 
K-hand luminosity of galaxies and, with a suitable choice of 
parameter values, it is possible to reproduce the bright end 
of the galaxy luminosity function. 



5 RESOLUTION EFFECTS 

The simulations discussed so far have a fixed mass resolu- 
tion, with the smallest resolved galaxies having a baryonic 
mass ~ W^'^h'^ Mq, roughly one tenth that of an L* galaxy. 



Resolution 


Feedback 


e 


T 


/gal 


^gal 


Low-res 


Thermal 


1.0 


1.0 


0.037 


11 


High-res 


Thermal 


1.0 


1.0 


0.053 


41 


Low-res 


Kinetic 


0.1 


N/A 


0.046 


11 


High-res 


Kinetic 


0.1 


N/A 


0.054 


36 



In the absence of feedback, resolution effects completely de- 
termine the outcome of a simulation: as smaller objects are 
resolved, increasingly large amounts of gas cool, raising both 
the global fraction of cooled material as well as the mass of 
the majority of the galaxies. Feedback limits the amount 
of gas that cools and so it seems plausible that an appro- 
priate feedback prescription might lead to results that do 
not depend strongly on resolution. In this case, an increase 
in resolution would lead to smaller galaxies being resolved 
without affecting the properties of the larger objects. 

To address the issue of resolution, we performed two 
additional simulations, using both thermal and kinetic feed- 
back models. For thermal feedback, we set e = 1 and r = 0.1; 
for kinetic feedback, we set e = 0.1. The cosmological pa- 
rameters and initial conditions were identical to those used 
for the SCDM model in except that the number of gas 
and dark matter particles were increased from 32768 (32^) 
of each to 110592 (48^). This results in a decrease in the 
minimum baryonic mass of a resolved galaxy by about a 
factor of 3, from ~ 1.6 x I0^°h~^ Mq to ~ 4.8 x lO^h'^ Mq. 
When constructing the initial density field, the phases of the 
large-scale waves were preserved, so that the same overall 
large-scale structure is present in all the simulations. We 
also increased the spatial resolution by decreasing the soft- 
ening length in proportion to m^^^, where m is the dark 
matter particle mass. This corresponds to a reduction by a 
factor of 1.5, for the duration of the simulation, to a final 
(Plummet) of ~ 6.7 kpc. All other parameters were fixed 
at the values used in the low resolution simulations. 

Table ^ compares the properties of the high resolution 
simulations with their low resolution counterparts. The last 
two columns give the fraction of baryons in resolved galax- 
ies, /gal, and the total number of galaxies, A'^gai. Galaxies are 
identified as before, using a friends-of-friends group-finder 
on the distributions of gas and star particles with a link- 
ing length of fe = 0.2, and only those with 32 or more star 
particles are retained. As expected, an increase in resolution 
leads to an increase in both the number of resolved galaxies 
and the galaxy baryon fraction. 

The main differences in the simulations with low and 
high resolution occur at early times and in small haloes. 
For example, the star formation histories (Fig. |ll| are quite 
similar except within the first gigayear, when smaller objects 
are resolved in the high resolution simulation. Increasing the 
resolution serves to extend the range of resolved galaxies 
to smaller masses, but has relatively small effects on more 
massive galaxies. For example, in Fig. ^ we plot the median 
ratio of cold gas to stars against stellar mass for galaxies 
with non-zero gas content and 32 or more particles (error 
bars represent 10 and 90 percentile points within each bin, 
of width AlogAf, — 0.5 for low resolution simulations and 
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Figure 11. Cumulative fraction of baryons in stars against time, 
in Gyr, for simulations with thermal and kinetic feedback, and 
different mass resolution. 



Figure 13. i^-band galaxy luminosity function for the simula- 
tions with thermal and kinetic feedback and varying mass resolu- 
tion. Also plotted is the best— fit Schechter function to the result 
from the 2dF and 2-Mass surveys (Cole et al. 2001). The value of 
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Figure 12. Median ratio of cold gas to stars against stellar mass 
for galaxies with non-zero gas content and 32 or more particles, 
in feedback simulations with varying mass resolution. Error bars 
represent 10 and 90 percentile points in each bin. 



0.35 for high resolution simulations). The cold gas fractions 
in the high resolution simulations are consistent with the 
low resolution simulations for galaxies resolved in both. 

The extension of the resolved galaxy population to 
smaller masses as the resolution is increased is clearly il- 
lustrated by the if-band galaxy luminosity functions, plot- 
ted in Fig. In the thermal feedback case, the luminosity 



function in the higher resolution simulation matches well 
on to the low resolution luminosity function. In the kinetic 
feedback case, there is some mismatch around A/^, but the 
number of galaxies is so small here that these differences 
are not significant. We conclude that, overall, the increase 
in the resolution of the simulation has little effect on the 
population of massive galaxies. 

Resolution affects are more noticeable in the amount 
and distribution of reheated gas. Fig. ^ shows the frac- 
tion of gas at different epochs lying in regions of different 
overdensity (c.f. Fig Increasing the resolution results in 
an increase in these fractions at all redshifts. The hydrogen 
column density distribution of reheated gas is illustrated in 
Fig.^ for the four simulations, at 2; = 3, 2, 1 & 0. Column 
densities were calculated by smoothing the density of local 
particles onto a 64 x 64 grid, using the SPH kernel with 
width set by the particles' smoothing lengths. Overlaid as 
circles (with radii proportional to M^!^^) are the positions of 
galaxies with 32 or more star particles. At z = 3 the gas in 
the low resolution simulations is confined to the few regions 
where the first galaxies have formed (most of these are unre- 
solved.) In the simulations with higher resolution the gas has 
a much larger covering factor because many more galaxies 
are resolved and so more of the volume is occupied by galax- 
ies. The range of column densities is high in these regions 
(reaching ~ lO^'^cm"^). At lower redshifts, the distribution 
of gas becomes much more skewed towards regions of higher 
overdensity, as the galaxies and haloes merge into larger 
units. By z = 0, the variation of density contrast in the gas is 
prominent, although the distribution is smoother in kinetic 
than in thermal models. As was inferred from Fig. kinetic 
feedback is able to transport gas into regions of lower den- 
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Figure 15. Column density of reheated hydrogen gas in SCDM simulations, with the positions of the galaxies superimposed. From top 
to bottom, the rows are for redshifts 3,2,1 & respectively. From left to right, the columns are for thermal feedback in the 32'^ & 48^ 

q q 1/3 

simulations, followed by kinetic feedback for 32'^ and 48 simulations. Radii of the circles are scaled proportional to M ' . 
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Figure 14. Fraction of baxyons supplied with supernova energy, 
at a given time, for thermal and kinetic feedback simulations with 
varying resolution (solid lines). Also plotted are subsets with 5 < 
500 (long-dashed lines), 5 < 50 (short-dashed lines) and <5 < 10 
(dotted lines). Redshift intervals are marked along the top of the 
figure. 

sity faster than thermal feedback, even though both cause a 
similar reduction of star formation by a = 0. 



6 CONCLUSIONS 

This paper has surveyed various phenomenological models 
for incorporating the effects of star formation and associated 
energy feedback (from supernovae) within cosmological sim- 
ulations of galaxy formation, using the SPH code HYDRA). 
Our main conclusions are as follows. 

(i) Most prescriptions for turning cold gas into stars 
used in previous SPH simulations give similar results when 
adjusted to give the same present day stellar density, n«. 
The mass in stars formed as a function of time is primarily 
determined by the gas density at which star formation is as- 
sumed to proceed efficiently. Thus, adopting a global density 
threshold (either in comoving or in physical coordinates), a 
Schmidt law, or a more complicated prescription involving, 
for example, the Jeans mass, makes little difference to the 
global properties of the galaxy population that forms. Given 
the limited resolution of current simulations (in which the 
internal properties of galaxies are unresolved), it is simplest 
to adopt a global density threshold. 

(ii) Incorporating star formation within a cosmologi- 
cal simulation makes no significant difference to the rate 
at which gas cools. The amount of gas cooled by a = 
only varies up to 1 per cent in simulations with and with- 
out star formation and the present day baryonic (gas+stars) 
mass function of galEixies also remains unchanged. However, 
when cold, dense gas particles are replaced by coUisionless 
particles, a saving of around 50 per cent of the CPU time 



is achieved from reducing the (computationally-expensive) 
work of calculating hydrodynamical forces. 

(iii) Simulations without feedback predict cooled gas 
fractions of around 35 per cent. Not only are these values 
resolution-dependent, but they exceed observational mea- 
surements by around a factor of 5. The problem of "over- 
cooling" is further highlighted in the iiT-band galaxy lumi- 
nosity function, which requires the majority of cooled gas to 
be "hidden" in low-mass non-visible stars (such as brown- 
dwarfs), in order to provide a reasonable match to the ob- 
served luminosity function. 

(iv) Injecting supernova energy to surrounding gas par- 
ticles, as thermal or kinetic energy, can reduce the star for- 
mation rate appreciably. However, thermal models can only 
achieve this if the reheated gas is prevented from cooling 
for a considerable fraction of a Hubble time. For maximal 
thermal feedback (i.e. all available energy is injected into 
the gas), the results are less sensitive to the choice of delay 
timescale than for models with a lower efficiency. Average 
fractions of cold gas in galaxies are unaffected by the inclu- 
sion of feedback. 

(v) A significant fraction of the gas reheated by super- 
novae is able to escape galaxy haloes and move into the low 
density regions {& < 10) responsible for the Ly-a forest. If 
this gas is viewed as a tracer of reprocessed material, this 
mechanism efficiently transports metals into low-density re- 
gions. 

(vi) The efficiency with which supernovae expel gas 
from haloes depends on the circular velocity of the halo 
( 100 — 300kms~^ in these simulations) - more gas is ex- 
pelled from smaller haloes, particularly for models capable 
of strong feedback. Larger simulations are required to estab- 
lish robustly the form of this relation. 

(vii) Increasing the mass resolution by about a factor of 
~ 3 increases the total amount of resolved galactic material 
primarily because smaller galaxies are then resolved. Conse- 
quently, the /C-band galaxy luminosity function is very sim- 
ilar over the range of magnitudes resolved in our high and 
low resolution simulations. The ratio of cold gas to stars is 
also essentially unaffected, as is distribution of gas perturbed 
by supernovae, although increasing the resolution increases 
the amount of this material, due to the higher number of 
resolved galaxies. Hence, until numerical convergence is es- 
tablished, the simulations underpredict the amount of this 
material. 

To fully tost the effectiveness of the feedback schemes 
that we have explored in this paper will require larger sim- 
ulations to be performed in statistically significant volumes. 
This will allow definite predictions to be made for e.g. the 
galaxy luminosity function and correlation function, which 
can readily be compared with observations. 
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